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Abstract. Using Beck and Cohen's superstatistics, we introduce in a systematic 
way a family of generalised Wishart-Laguerre ensembles of random matrices with 
Dyson index (3 = 1,2, and 4. The entries of the data matrix are Gaussian random 
variables whose variances 77 fluctuate from one sample to another according to 
a certain probability density f(ri) and a single deformation parameter 7. Three 
superstatistical classes for f{r]) are usually considered: x^-i inverse x^- and log-normal- 
distributions. While the first class, already considered by two of the authors, leads 
to a power-law decay of the spectral density, we here introduce and solve exactly 
a superposition of Wishart-Laguerre ensembles with inverse ^^-distribution. The 
corresponding macroscopic spectral density is given by a 7-deformation of the semi- 
circle and Marcenko-Pastur laws, on a non-compact support with exponential tails. 
After discussing in detail the validity of Wigner's surmise in the Wishart-Laguerre 
class, we introduce a generalised 7-dependent surmise with stretched-exponential tails, 
which well approximates the individual level spacing distribution in the bulk. The 
analytical results are in excellent agreement with numerical simulations. To illustrate 
our findings we compare the x^- and inverse x^-class to empirical data from financial 
covariance matrices. 



1. Introduction 

Random matrix theory (RMT) is known to find applications in many physical systems 
[H [2] . Its central assumption is that the Hamiltonian of the system under consideration 
can be replaced by an ensemble of random matrices that are consistent with its 
global symmetries. The matrix elements of such a random Hamiltonian are typically 
independent Gaussian variables with mean zero and variance one in the real, complex or 
quaternion domain. These are respectively labelled by the Dyson index /3 = 1, 2 and 4, 
which in turn corresponds to the invariance group of the ensemble (orthogonal, unitary 
or symplectic). Another set of ensembles of random matrices, which is well studied in 
RMT is known as Wishart-Laguerre or chiral ensemble and will be denoted by WL in 
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the following. The WL ensemble contains random matrices of the form [3] W = X^X, 
where X is a rectangular matrix of size M x N [M > N), whose entries are independent 
Gaussian variables in the simplest case, and X''' is the Hermitian conjugate of X. As 
such, the WL ensemble contains (positive definite) covariance matrices W of maximally 
random data sets. They have since appeared in many different contexts ranging from 
mathematical statistics, statistical physics and quantitative finance to gauge theories, 
quantum gravity and telecommunications [1]. 

The Gaussian distribution of matrix elements above can be obtained by maximising 
the Boltzmann-Gibbs-Shannon entropy subject to the constraint of normalisation and 
of fixed expectation value of Tr(X"''X) This observation has been taken by 

several authors [6] as a starting point for generalising the classical Wigner-Dyson (WD) 
ensembles, by maximising Tsallis' non-extensive entropy subject to the same constraints. 
The resulting distributions of matrix elements are no longer Gaussian but rather follow 
a power law (for a different generating mechanism of such ensembles see [7]). Recently, 
a similar generalisation of the WL ensembles has been worked out [8], following earlier 
studies [9] on the covariance matrices of financial data. 

In this context, it is known that the spectral density of covariance matrices 
built from empirical data may deviate significantly from their purely random (WL) 
counterpart: in order to obtain a better agreement, it is necessary to introduce 
correlations among the matrix elements of X, typically allowing the variance of each 
entry to fluctuate. In [10] this is done using a multivariate Student distribution, where 
the random entries of X are written as a product of two variables with a Gaussian 
and inverse x^-distribution respectively. Generically, this approach spoils the complete 
solvability of the model and prevents from going beyond the average spectral density 
[To] . In order to study other correlations, it is therefore of great interest to introduce 
generalisations of WL where the independence of X-entries is dropped, but the exact 
solvability is retained. One of the first examples of such models was provided in [8], 
where a good fit to the power-law decay of financial covariance spectra was obtained. 

It is desirable to place the model [8] within a more general framework, as was 
done previously for the WD class [Til [12] • The key observation is to resort to the 
ideas of superstatistics (or statistics of a statistics) proposed by Beck and Cohen [13] . 
Outside RMT, this formalism has been elaborated and applied successfully to a wide 
variety of physical problems, e.g., in [T3]. In thermostatics, superstatistics arises as a 
weighted average of ordinary Boltzmann statistics due to fluctuations of one or more 
intensive parameters (e.g. the inverse temperature). Typically, the distribution of the 
superstatistical parameter falls into three 'universahty' classes: x^-y inverse x^- ^-nd log- 
normal distributions (see section [2TT] for a more detailed discussion). Superstatistical 
RMT [11] analogously assumes that the Hamiltonian of the system is locally described 
by a standard WD ensemble with a given variance, and when averaging over the 
whole system the variance is integrated over with a specific distribution. Consequently, 
superstatistical RMT is a superposition or integral transform of the usual Gaussian 
RMTs. 
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The same concept of integral transforms of standard RMT has also appeared in 
other contexts, e.g. the fixed or restricted trace ensembles [15], also called norm- 
dependent ensembles. Very general results can be deduced for them, without specifying 
the distribution [121 IIS]- However, no systematic applications of superstatistics and 
integral transforms to the WL class are known so far. It is the purpose of this paper to 
fill this gap. 

We are going to introduce first a family of generalised WL ensemble depending on 
a single real parameter 7, where the prescriptions of superstatistics are incorporated 
in a systematic way, and then solve exactly a representative example where the entries 
of X are Gaussian random variables whose variances fluctuate from sample-to-sample 
according to an inverse x^-distribution. Given the appearance of such a distribution in 
modelling the volatility of financial markets [10], it is sensible to compute the tail of the 
spectral density in our model and compare it to the findings in [10] and [8]. However, 
because of WL ensembles being tailored to times series of data, we expect applications 
beyond financial covariance matrices as applied here. 

Despite the complicated correlations among the entries, and thanks to the integral 
transform mapping we can go to an eigenvalue basis and write the joint probability 
density function (jpdf) of eigenvalues as an integral over the corresponding WL one, a 
feature that was already exploited in [8] for the x^-distribution. Being able to perform 
all integrals analytically, we have full control over exchanging the large A^-limit with our 
deformation parameter 7. In the limit 7 —>■ 00, our generalised ensembles are designed 
to recover the standard WL, and this fact will constitute an important consistency check 
in the following. A third class of models can be obtained by folding with the log-normal 
distribution, but we will not deal with this case here, lacking a full analytic solution. 

It is worth mentioning that the applicability of this formalism is not limited to the 
above three universality classes. In particular, it would be very interesting to investigate 
whether the choice of the best superstatistical distribution could be inferred from the 
data set, instead of being somehow 'postulated' a priori. This would be reminiscent of 
the Bayesian approach to superstatistics [IT], this time in the more complicated RMT 
setting. 

In the next section [2], we briefly review the prescriptions of superstatistical RMT, 
together with the three universality classes in subsection 12. 1[ Then we give a step- 
by-step derivation of our superstatistical WL ensemble in subsection 12.21 In section 
E] we define the spectral properties to be computed in our model and discuss their 
universality. We then investigate the large matrix size limit of the spectral density 
in section HI where we distinguish between square and rectangular matrices of size 
= cM in the two subsections. In section [5] we first present a detailed discussion 
on the applicability of the standard Wigner's surmise for the level spacing in the WL 
class, and then we derive a new, generalised 7-dependent surmise which is appropriate 
to describe with excellent approximation the individual level spacing in the bulk of 
our superstatistical model. Numerical checks on the analytical results are provided 
throughout the text and the algorithm used is described in [Appendix A[ In section [6] 
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we compare two superstatistical distributions to empirical data from financial covariance 
matrices, before offering concluding remarks in section [71 

2. Superstatistical Wishart-Laguerre ensembles 

The probability density of the data-matrix entries in WL ensembles is given by 

PwLiv,^) = ^^exp [-r//3Tr(XtX)] , (1) 

where rj is proportional to the inverse variance of matrix elements. 
The normalisation is given by the partition function 

ZMV) = y «exp [-„T.(X.X)] -^{-^) . (2) 

As mentioned above, X is a matrix of size M x N with real, complex or quaternion real 
elements for the values /5 = 1, 2 or 4, respectively. We parametrise N = cM for later 
convenience, where c < 1 distinguishes two different large-A^ limits. The integration 
measure dX is defined by integrating over all independent matrix elements of X with 
a flat measure. The statistical information about the positive definite eigenvalues of 
the Wishart matrix W = X"''X, or equivalently the singular values of the matrix X 
can be obtained integrating out all the undesired variables from the distribution of the 
matrix elements, using its orthogonal, unitary or symplectic invariance for /3 = 1, 2 and 
4 respectively. 

We now define a superstatistical family of WL ensembles. The key ingredients are 
the following: 

(i) A real deformation parameter 7 > 0, which roughly quantifies how far the new 
model lies from the traditional, unperturbed WL ensemble. In the limit 7 — > 00, 
we expect to recover WL exactly. 

(ii) A normalised probability density /{f]), such that 

POO 

1= / dvfiv) . (3) 
Jo 

It is understood that /(t^) depends on 7 as well, but we will not show this 
dependence explicit in order to keep the notation light. 

The probability distribution of matrix elements for the generalised model is then 
obtained as follows: 

P(X) = drj f{v)PwL X) = {PwL {iilH X) (4) 

where (■)/ means average over the distribution /(?7) and £(7) is a simple function of the 
deformation parameter (see subsection 12.21 below for details). The choice of distribution 
f{ri) is determined by the system under consideration. In the absence of fluctuations of 
the variance, f{r]) = 5{rj — 770) and we reobtain the standard WL ensemble. Typically, 
the distribution firf) depends explicitly on the parameter 7, in such a way that for 
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7 oo a delta- function limit is obtained and thus the WL results (before or after 
taking large) are duly recovered. 

While the distribution defined in eq. (HI) is formally normalised, when exchanging 
the integration J c/X with J drj, the prescription given so far is not complete. The choice 
of a distribution /(?]), and in particular of the iV-dependence of its parameter (s) has 
to be such that a) the integral over -P(X) is convergent, and that b) an N independent 
limit for the spectral density can be found after a proper rescaling of variables. We will 
come back to this issue below. 



2.1. The three superstatistical classes 



Beck et al. [18] have argued that experimental data can be described by one of three 
superstatistical universal classes, namely the x^-? inverse x^-? or log-normal-distribution. 
Below we briefly discuss each class and the properties of its corresponding distribution, 
before turning to the detailed solution for the inverse x^-class in the next sections. 



1) X 



^ distribution: 



The x^-distribution with degree u and average rjo is given by 



r(|) 



2^0 



rj'- 



exp 



— urj 



2vo 



(5) 



where we define rjQ = rif{ri)dri = {rf). This distribution is appropriate if the 
variable t] can be represented as a sum of squares of v Gaussian random variables. 
The superstatistical distribution arising from ^ is Tsallis' statistics with power 
law tails [19] and is believed to be relevant e.g. for cosmic ray statistics [20] . 

2) inverse x^-distribution: 

This distribution is found if ?7~^, rather than 77, is the sum of several squared 



Gaussian random variables. Its distribution /2(?7) is given by 



r(i) 1 2 ; 



exp 



(6) 



with degree u and average 779. 

This class is appropriate for systems exhibiting velocity distribution with 
exponential tails [21] and was shown to describe the spectral fluctuations of billiards 
with mixed regular-chaotic dynamics better than the other two distributions |22j . 
3) log —normal distribution: 

In the third class, instead of being the sum of contributions, the random variable 
r] may be generated by a multiplicative random processes. Then its logarithm 
log(?7) = ^^=ilog(xj) is a sum of z/ Gaussian random variables Xj. Thus it is 
log-normally distributed: 



fsiv) = 
It has an average (77) 



1 



2-7r VT] 



exp 



2v^ 



(7) 



fI^/w and variance fi'^w{w — 1), where w = exp(t;^ 
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This class has been found relevant for Lagrangian and Eulerian turbulence 
Because apparently this class does not lead to closed analytical results for the 
distribution of matrix elements -P3(X) or its correlation functions we do not use it 
as an example in this paper. 

2.2. Building the superstatistical ensembles 

Since the prescription (jlj) is not completely straightforward, we provide here a step-by- 
step construction of the superstatistical probabihty density -P(X). We focus here on 
the inverse x^-distribution (I6j) as a new example. 

(i) We start with the usual WL probability density for the entries, eq. ([1]): 

P (X) exp [-r//3Tr(XtX)] . (8) 

Henceforth we omit the normalisation constants. 

(ii) Next, we convolve the WL weight with the inverse x^-distribution iQ: 



exp [-r//3Tr(XtX)] . (9) 



(iii) The u parameter becomes the deformation parameter of the model as 7 = z//2. For 
a clearer notation, we now make the further replacement rj = and eventually 
integrate over the possible values for 

P(X) - ^"^e^exp (^--^ exp [-e7/?Tr(XtX)] . (10) 

Note that the average ?7o is no longer needed explicitly and has been absorbed in 
the other parameters. In fllOp . the deformation parameter 7 appears in both the 
superstatistical weight function and inside the WL Gaussian weight (here ^(7) = 7, 
while in the [S] model £(7) = I/7). 

(iv) Although formally exact, eq. ffTOj) does not lead to a A^-independent spectral 
density, as can be quickly realized with a modest amount of foresight (see 0211) 
and (122|) ). The last step is thus to amend slightly (fTOj) . including a suitable extra 
factor for normalisation ^-(/3/2)A'Af. 

(X) ^ ^^P ^^P [-e7/3Tr(XtX)] . (11) 



(v) Eventually, the integral in ffTTj) can be evaluated, and gives the probabihty density 
for the entries of our superstatistical WL model as: 

/ , N i (7+1-/3 ArM/2) / N 

P2 (X) cx (Tr(XtX)j K^+,_pMM/2 (2vWlMXtX)j . (12) 

Here K^{z) is the modified A'-Bessel function of second kind. 

This distribution is the main subject of this paper. It is convergent for all values of 
7 > 0, and we will show later that our choice of the normalisation ^-(/5/2)A^Af effectively 
leads to an A^-independent limit for the average spectral density at fixed and finite 7. 
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We can analytically recover the standard WL ensemble for 7 — oo (at fixed and 
finite A^, M) using the following non-standard asymptotic of the Bessel-ZT functior(||: 



^2^^ '^exp[ 



(13) 



While we did not find eq. (11311 in tables, it follows easily form the following integral 
representation 8.432.6 pi], after using a saddle point approximation 



X ' 



27+1 



dt t-^-^ 



exp 



X 

-'-li 



x' >0 



(14) 



To derive eq. ( |T3i) we have included the fluctuations around the saddle point and used 
the standard notation / ~ to mean that f /g — > 1. The 7-dependent prefactor can be 
cancelled by a proper normalisation. 

It is important to stress that the choice of the normalisation ^-(/3/2)MAr ^Qgg 
affect the correct 7 — > 00 asymptotics in any way, while being the only sensible 
prescription when taking the large A^, M limit at fixed 7 for the average spectral density 
(see section Hj). 

It is also worth mentioning that, following the same steps as above, one can 
incorporate the other two superstatistical distributions. For example, for the x^- 
distribution, one is naturally led to: 



Pi(X) oc 



d^ e-i+^^^^^e-«exp 



- ^ Tr(XtX) 

7 



oc 



7 



TrfX+X^ 



--f+^f3NM 



(15) 



The power-law decay of the matrix elements translates into all correlation functions [8], 
depending on a single rescaled parameter, = 7 — ^NM — 1 > 0, which is kept fixed in 
the large- iV limit. The reduction back to standard WL can be made both before (and 
after) the large- limit using 



lim (1 + 7 

7^00 



(16) 



For a more detailed discussion we refer to El. 



3. Generalisation of WL with exponential tails 



After introducing the probability distribution of matrix elements in our generalised 
inverse x^-WL eq. f fTTj) we go to an eigenvalue basis in this section and define all 
correlation functions. 

The corresponding joint probability distributions function (jpdf) of positive definite 
eigenvalues Ai, Aat of the matrix W = X"''X reads 



oc 



d^ 



exp 



(17) 



I Note that the argument y^z of the Bessel-A' function becomes quadratic in the exponent. 
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Kj+i-l3NM/2 



A 




Here we have used the jpdf Vwl of the standard WL ensemble, which is given by 

N N 

P^,(A,,...,A^;OocniA,-A.|^n4^*™)-V«-^^-\ (18) 

i<j k=l 

For convergence 7 has to be positive. Both jpdf's still have to be normalised by the 
respective partition function and Zwl{,C)^ obtained by integrating over all eigenvalues. 

The fc-point density correlation functions defined by integrating out N — k 
arguments of the jpdf are given as: 

-R7,fc(Ai, Ayfc) = — J d\k+i---J d\NVj{\i, \n)- (19) 

Because of the linear relationship between the jpdf of our ensemble and WL, eq. ( flTl) . 
we can express the fc-point functions of the former through the latter: 

-R7,fc(Ai, Afc) = [ dC, — , „ p.^.. e g '^^^'^ -RvKL,fc(Ai, Afc; ^), (20) 
Jo ^1+2- ^NM 

where -RiyL,fc(Ai, A^; (^) is the corresponding /c-point function for the WL ensemble 
defined in the same fashion. The ^-dependent ratio of the two partition functions which 
are also linearly related through 

/.oo 2 

^7 = / ZwdO , (21) 



^7+2- 



easily follows from eq. ([2]): 



NM 



~ r(7 + i) • ^^^^ 

Comparing eq. ( l22i) and ( 120|) , it is already apparent that the choice of the normalisation 
in eq. (fTOj) was in fact necessary to neutralise exactly the A^- and M-dependence coming 
from the partition functions. In all our formulas, the limit 7 — > 00 at finite TV and M 
leads back to WL, using eq. (fT3ll . 



3.1. Universality 

We close this section by discussing the universality of generalisations of WL. Following 
from a variational, generalised entropy principle we should expect a certain robustness 
of our results, which is indeed the case. 

First, all fc-point correlation functions of the WL ensembles are explicitly known for 
finite and infinite- iV, in terms of Laguerre polynomials and their asymptotics. Thanks 
to the integral mapping eq. fl20|) . the superstatistical models are exactly solvable as well 
(see e.g. [12], [8]). Also because of this, we could allow for a more general WL ensemble 
with a polynomial potential V instead of the Gaussian, Tr(X"''X) Trl^(X''"X), in 
order to probe the universality of our results under deformations. These non-Gaussian 
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ensembles can again be solved using the technique of (skew-) orthogonal polynomials 
for any finite N and M, and again the complete solvability of the superstatistical 
generalisation is guaranteed by the very same integral relation eq. ( l20l) . 

It is known that in the microscopic large-A^ limit the asymptotic of the 
orthogonal polynomials for these non-Gaussian models is universal [251126] (for rigorous 
mathematical proofs see [27] and references therein). This implies that the microscopic 
correlations in our model remain unchanged after the integral transform. As an example 
for such a microscopic quantity we deal with the level spacing distribution in section 
|5l For a detailed discussion of the universal microscopic spectral correlations with a 
X^-distribution we refer to [8], including references. 

About the macroscopic large-iV limit, the (generalised) semi-circle and Marcenko- 
Pastur densities are known to be less robust. They are altered by deformations via a 
polynomial potential V but still remain computable, see e.g. the discussion in [15] and 
references therein. However, a different kind of deformation is known that leaves the 
semi-circle or Marcenko-Pastur density unchanged: the so-called Wigner ensembles (see 
[28] for a review). The Gaussian random variables of the WL ensemble are replaced 
by independent random variables with zero mean and fixed second moment. While 
this generalisation destroys both the invariance as well as the integrability of higher 
correlation functions, its large- macroscopic density is the same as in WL, and thus 
also our integral transform of WL, after taking the large- A^ limit. 

4. The macroscopic spectral density 

In this section we focus on the simplest observable, the spectral density R^^k=iW 
obtained by integrating out all eigenvalues but on^. 

When taking the large- A^, M limit, we keep the combination c = N/M < 1 fixed. 
This leads to two different behaviours for the WL spectral density: the semi-circle (in 
squared variables) for c = 1 and the Marcenko-Pastur (MP) density for c < 1, see eqs. 
( !24|) . (138|) below. Therefore we expect two different limits for our generalised ensembles 
as well, and the generalisations of these two well-known WL results will be dealt with 
in two separate subsections below. 

Although the average spectral density for WL (and hence for our model as 
well through eq. ( l20i) ) is known explicitly for any finite- A^, M in terms of Laguerre 
polynomials (see e.g. [21 [8]), the generalised spectral density for large A^ can be obtained 
following a much simpler route: we replace the finite- A^ WL quantity under the integral 
by its large- A^ result. The correctness of this approach has already been shown in [8]. 
As a further check, we can eventually take 7 — * oo to recover the semi-circle or MP 
density. 

§ Following the general definition cq. ([TO]). i?(A) denotes the spectral density normalised to N. 
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4.1. Generalised semi-circle for law c = 1 

In the case c = 1, the large- asymptotic expression for the density of eigenvalues of a 
WL ensembles with distribution exp[— 77/? Tr X"''X] eq. ([I]) is given by 



Ion 

lim RwL i(A) = -W^ - 1 with A G (0,2Ar/r/] , (23) 
Af>i ' 71 V rjX 

and otherwise. An A^-independent macroscopic density with average one is obtained 
by rescaling A — > x {X)wl with its mean {X)wl = Tr X"l'X)vt/j^ = calculated with 
respect to the above distribution. We obtain 



pwLix) = lim l(A)i?H/L,i(x(A)) = ^J--l, with xe(0,4].(24) 

It diverges as \l\fx at the origin and vanishes as a square root at the upper edge of 
support. Eq. (12^ corresponds to a semi-circle after mapping eigenvalues from M+ to 
M, and we will comment more on that below (see Fig. [1]). 

Following the same procedure as above for the inverse ensemble, we first need 
to determine the average eigenvalue in our ensemble (see e.g. [8] appendix A): 

1 /"^ 1 1 M('7-hl) 

^T,l '^^;^rp7n^-' 2-40{a«)> = , (25) 

using eq. fl22l) and the mean in WL above. 

Next we define the generalised macroscopic A^-independent density as 

p^(x) = lim ^{X)^R^{x{X)^) 

N—*oo iV 



. 1 (A), f,^ \ _ 1 . (26) 

Here we have simply inserted the large- A^ WL density eq. ( 1231) with the parameter 
?7 = ^7 into the integrand eq. fl20|) . Consequently the interval of integration is truncated 
and given by X = (0, ^.^^^'^j^^ ]. Thanks to the correct choice of normalisation in eq. ffTTj) . 
the A^-dependence completely drops out and for any fixed 7 > the new spectral density 
admits the following integral representation, after changing variables t = — 



, , (7 + 1)^+^ (x\ . , , 
"^^"^ = 2vrr(7 + l) U) / ^'^"P 



{t + iy~^Vt. (27) 



4 

This is our first main result of this section. This density as well as its first moment are 
correctly normalised to 1: 







p^{x)dx = 1 = xp^{x)dx . (28) 
'0 Jo 

As a consistency check, we can now take the limit 7 ^ 00. The integral becomes 
amenable to a saddle point approximation, and taking into account the fluctuations 
around the saddle point we reproduce exactly eq. ( l24l) as we should. Eq. ( 127|) is plotted 
in Fig. [1] (left) for various values of 7. Our generalised density has support on the full 
M+, in contrast with the compact support of the WL semi-circle. 
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In order to analytically derive the asymptotic behaviour of our generalised semi- 
circle eq. (1271) at the origin x and at infinity x ^ oo it is more convenient to 
express the integral through special functions. Using eq. 3.383.5 [2^ we can express the 
density eq. (127|) as 



(7 + 1)^"^^ 



exp 



-(7 



4) * 2-^' 



3 / \ 

2; (7 + 1)4 



(29) 



4v/^r(7 + l) 

where \&(a,6;z) is the Tricomi confiuent hypergeometric function (sometimes denoted 
U{a, b] z)). The asymptotics for x — > is easy to obtain as: 

(7 + 1)^/^1(7 + 1/2) 1 



X 



(30) 



vrr(7 + 1) 

A few comments are in order. First, the 7-dependent prefactor in (!30l) reproduces 
the correct WL asymptotic behaviour from (l2ll) Pwl{x) ~ for x — > when 7 — s> 00. 
Secondly, the WL inverse square root divergence at the origin is not modified in the 
superstatistical generalisations, a feature that is shared by the x^-deformed WL [8]. 




Figure 1. The generalised semi-circle eq. ((29)) before (left) and after (right) the 
mapping eq. ((33)) to K for different values of 7 = 0.5, 2, 10 (green, blue, violet 
respectively). For comparison, the standard WL semicircle is plotted as well (black). 



Next we turn to the large-x asymptotics. Using only the leading order term [30] 

-^{a,b;z) ~ z'"" for j^] -> 00 , (31) 
we obtain immediately the following result from eq. (129)1 

^ 40Fr(7 + i) ^"p r^^ + (4) • ('2) 

We illustrate the behaviour of the generalised semi-circle density eq. (1291) in Fig. [T] In 
order to better resolve the behaviour at the origin, we map our spectral density from 
the positive real axis to the full real axis by defining 

= \y\ P,{y^) , (33) 



Superstatistical Wishart-Laguerre Ensembles 



12 



that is to a normalised density on M, dy'd^{y) = 1. Furthermore, this map leads the 
second moment to be normalised dyy'^'d.yiy) = 1. In particular we obtain for WL a 
semi-circle from eq. (l24l) . 

^WL{y) = ^ on [-2, 2] . (34) 

In this squared- variables picture, a further, interesting feature of the generalised spectral 
density eq. (!27|) appears when considering the opposing limit 7^0. There, the spectral 
density admits the following simplified form: 

r 



2 ' 4 

d^^Q{x) = \x\ \ ^ (35) 

(where r(x, y) is an incomplete Gamma function), which displays a cusp at x = 0. This 
feature often appears in spectral densities of complex networks |31j . 

As a final remark in the standard WL (or WD) class the density also has an 
exponential tail at finite- iV, exp[— A^x^]. However, in the macroscopic large- A^ limit 
it disappears while in our deformed model such a tail persists. 

4-2. Generalised Marcenko-Pastur law for c < 1 



In this subsection we deal with the limit in which the matrix X remains rectangular, 
that is both M and A^ become large with N/M = c < 1 fixed. We follow the same 
steps as in the previous subsection, first recalling the results for WL that need to be 
incorporated into our model. It is known that in the standard WL ensemble eq. ([I]), 
the average density of eigenvalues is given for large- A^ as follows: 




lim RwlAX) = { A - i-X^ ) ( - A ) , with A G 



'^x ^x' 



.(36) 



Af>l TtA 

Here we have defined the following bounds for later use 

= {c-i ± If , with < c < 1 . (37) 

In the limit c — > 1 we recover from eq. (l36l) the semi-circle eq. (123!) from the last section. 
An A^-independent density with mean one is again obtained by rescaling with the mean 
eigenvalue position {\)wl = and normalising 

Pmp{x) = lim ^{X)wlRwl,i{^Wwl) = 7r~Vi^ - cX_)(cX+ - x) , (38) 
N-^oo i\ Zncx 

with X G [cX_,cX+] . 

This is called Marcenko-Pastur (MP) law [29]. Taking c — > 1 we again recover the 
A^-independent semicircle eq. ( 1241) . 

Turning to our generalised model we now insert the large- A^ result eq. (!36|) into eq. 
( !20l) and rescale with the mean eq. ( |25l) 

T^W^^^i^Wl) (39) 
N,M—^oo iV 
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where X _ , / , 

x(7+l) ' a:(7+l) 

eq. (122]) and substituting ^ 



is the truncated integration range. Inserting the ratio 
= ^x(7 + l)/c we arrive at the following 



X ' 



X] 



27rr(7 + 1) V c 



7 



7+1 



X4 



X- 



dt 



exp 



x(7 + 1) 
tc 



v/(t-X_)(X+-t) , (40) 



the second main result of this section. It is valid for any fixed 7 > with < c < 1. 
As a check we can take the limit 7 ^ cxd, and a saddle point evaluation including the 
fluctuations around the point = f leads back to the MP density eq. ( l38i) . 

Next we turn to the asymptotic analysis for x and x 00. For small values 
of X we simply obtain 

p^(x) ~ D^x'^ (41) 

where the constant easily follows from eq. fHOj) . For x — > 00, we obtain: 



^ 4v/^r(7 + i)(x+)7-i V c ; ^ ^ 

As a check, we can recover the correct asymptotical behaviour ( l32i) when c — > 1 (in this 
case X_ — > and X+ — > 4). Note that the decay for large arguments is the same as for 
the generalised semi-circle eq. ( I32l) . in terms of variables x/{cX^). An example for eq. 
fHOj) is shown in Fig. [2l Such exponentially decaying correlations occur in exponentially 
growing complex networks, apart from more frequent power law decays |31j . 



0/0 (^+1) - 




Figure 2. The generalised MP density eq. (|40l) for c 0.4 and different values 
of 7 = 0.5,2,10 (green, blue, violet respectively). For comparison, the standard MP 
distribution eq. p8|) is plotted as well (black). 



In addition to our analytical checks comparing to WL and c = 1 we have also performed 
numerical simulations, generating matrices in our generalised class with the algorithm 
given in Appendix A In Fig. [3] we compare the numerical results for the spectral 
density at finite N, M with the theoretical prediction (valid at infinite TV) eq. (140|) . We 
find an excellent agreement already for moderate N and M [N = 10, M = 40). 
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Figure 3. Numerical clicck for /3 = 1: Histogram of eigenvalues for N = 10, M ~ 40 
and 7 = 1, 10 (blue and magenta, respectively), together with the theoretical result eq. 
(|40ll for N = oo (black and red solid lines). The average is obtained over R = 50000 
samples. 



5. Level Spacing from a Wigner Surmise 

The exact computation of the level spacing distribution for finite matrix size (and in 
the microscopic large-iV limit in the bulk) is mathematically quite nontrivial, even in 
the simplest case of the WD ensembles, see e.g. [2]. 

Therefore long ago Wigner came up with the idea to compute the exact spacing 
distribution in the WD class for the simple 2x2 case, obtaining 



The /3-dependent constants a^, hp easily follow by fixing the norm and first moment to 
unity (see e.g. in [1]). Below we will mostly need (3 = 1 with oi = 7r/2, and hi = n/A. 
The result eq. (H3|) turns out to be an excellent approximation of the true result for 
moderately large A^, see e.g. Fig. 1.5 in [2], and we aim at a similar approximate 
representation for the spacing distributions in superstatistical ensembles. 

However, before undertaking this more complicated task, we ask the simpler 
question: does a Wigner surmise using N = 2 work well in the unperturbed WL 
ensemble? We have not seen this issue discussed in depth in the literature (see however 
[32] section V and references therein) and we feel it is appropriate to spend a few words 
on it in the next section. 

5.1. A Wigner surmise in WL? 

Starting from eq. ( fTSl) for N = 2 with weight exp[—nf3X\ [n = 1/2 for standard normal 
entries), we can compute the exact spacing distribution Vp{s) in WL by integrating over 
both eigenvalues with the constraint 6{X2 — Xi — s). Since this was already done in [8] 




(43) 
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we can just quote the result: 

)(.) = Cs^^'^'/'K^/,^,{nPs) , (44) 

where K^{x) is a modified Bessel function, i/ = f (M — TV + 1) — 1 = f (M — 1) — 1, and 
the constant C is given by: 

C = (^2-V2+/3+'^(n/3)-3/2-/3-'^r (^-^) ^ + + f ) ) ' ^^^^ 

to ensure a normahsation to unity. The first moment can be normahsed by computing 

d= dssVf\s) , (46) 
Jo 

and then defining 

Vf{s)^dvi^\sd) . (47) 

The result eq. P7|) is obviously different from eq. fj43l) . even for M = N = 2. Can 
eq. (1471) be a better approximation than eq. ( l43l) for the true WL spacing at finite 
but large? 

In order to compare, let us pick N = M and /3 = 1 or 2 implying u = or 
respectively. The resulting Bessel--ft' function simplifies only for half-integer index, 
-^1/2(2^) = \/7r/2xe~^, and we obtain 

Vl^_='l{s) sKois) , (48) 
2 

Viit\s) r-. exp[-2s] . (49) 

After normalising the first moment, neither case matches with eq. (143|) and for N ^ M 
(or /3 = 4) the difference is even more pronounced. 

In order to check which spacing distribution is a better approximation for larger 
instances, we performed numerical simulations using the Dumitriu-Edelman tridiagonal 
algorithm [33] for the case /3 = 1 and A^ = M, with the result shown in Fig. |H 
Before commenting on this we mention in passing our procedure. Usually, the numerical 
quantity to be compared to the surmise is the so-called nearest-neighbour spacing 
distribution, i.e. the normalized histogram of all the spacings among consecutive 
eigenvalues in the bulk after unfolding. This procedure is often time-consuming and 
special care is needed to avoid spurious effects in the analysis (see [321 [M] and references 
therein for a detailed discussion on unfolding procedures). 

We can thus resort to the following, equivalent method to extract the individual 
spacing at a given location k in the spectrum (see [35]) as: 

The average (■) is taken over many samples and obviously (sfc) = 1. For a given A^ we 
have picked different values of k to check that our results do not depend on the position 
in the bulk. 




Figure 4. Numerical clieck of tlie WL and WD surmise for /3 = 1 and = M: 
N ^ 2 vs. = 10 (13), for various locations in the bulk: k ~ 5, 7, 8. 



While the numerical histogram for N = 2 follows the exact result for the spacing from 
eqs. (l48l) . for increasing it quickly converges towards the WD surmise eq. (H3l) . We 
have checked that this holds when choosing N ^ M a.s well. 

What is the explanation? The reason lies in the well known fact that in the bulk 
both in the WD and WL classes the correlations are governed by the Sine kernel (for a 
rigorous proof, see [27]). The exact spacing distribution can then be expressed in terms of 
Fredholm eigenvalues of this kernel, see. e.g. [2], which happens to be well approximated 
by the WD surmise. This is why the WD surmise applies to both ensembles. In contrast 
the N = 2 WL surmise is not a good approximation, and so the surmise does not work 
here at all. In the next subsection, we seek for a generalised Wigner's surmise holding 
in the bulk of our superstatistical model. 



5.2. Level spacing for our generalised model 

How should an approximate spacing distribution for the generalised WL model with 
an inverse x^-distribution be constructed? Because of the failure of an = 2 surmise 
in WL (as we just explained), we should not expect that inserting it into our integral 
transform will work in the superstatistical case. Instead, we will follow a more heuristic 
approach confirmed by various numerical checks. 

Rather that following the procedure described in subsect. I3.1l we directly start from 
an integral transform of eq. (H3l) but with a .^-dependent varianc^ 

PSk^;0 = 2e^+^(|)'^r(^)"Vexp[-/3eV/2] . (51) 

II The quadratic rather than linear power in ^ in the exponent can be motivated by a change of variables 
in eq. (HH]) e~^'^ ^ from the WL to the Gaussian WD weight. 
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Eq. (1511) is normalised to one (the first moment will be normalised below). Because of 
/ d^fiO — 1 ^he following folded surmise for our generalised WL is also normalised: 



with 



d^r'^'e-lFfj,{s;0 = Cy 



1\-1 



exp 



(52) 



(53) 




Figure 5. Level spacing distribution eq. (|55p at /3 = 1,2,4 (from left to right) and 
different values of 7 = 1, 10, 50 (red, blue, violet respectively). We have also included 
the WD surmise eq. (j43p in green for comparison. 



To normalise the first moment we compute 

(7 + i)r(f + 1 



/3\ ^pZ/g+i 



(54) 



2 J \ 2 

and obtain the approximate spacing distribution in our generalised WL with an inverse 
X^-distribution: 

P!^^\s) = d^Pi^\sd,) . (55) 

Eq. fl55|l is the main result of this subsection. Given the known universality of the 
(approximate) eq. ( I43l) . our new spacing distribution will also be universal under a large 
class of deformations as discussed in subsection 13. 1[ The integral can be evaluated in 
terms of hypergeometric functions, but we prefer to keep the integral form for simplicity. 
We have checked explicitly that in the limit 7 — > 00 we correctly reproduce eq. (H3l) 
valid for WL (A^ > 2). 

Eq. ( I55l) is plotted for all three values of [3 = 1,2,4 and various 7 in Fig. [5l 
including the limiting WD distribution eq. (143|) . 

Our new surmise eq. (l55l) has the following asymptotic behaviour. The level 
repulsion at short distances s ^ is given by 

Kis™"^^'^) , (56) 

At large distances s ^ 00 however we get a new behaviour in terms of a stretched 
exponential 



(57) 
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which can be obtained after changing variables and making a saddle point approximation 
{a J = In both formulae, we have omitted the exact form of the 7-dependent 

prefactors ki,2- 
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Figure 6. Comparison to numerical simulations; top: for fixed 7 = 7 a comparison 
between eq. (|55p and numerically generated spacings for various combinations of N, M 
and locations k within the bulk; bottom: for fixed TV = 10, M = 15 eqs. ((55|) . (|43|) 
and numerically generated spacings for various values of 7 and locations k. 



To test our new surmise, we have performed numerical simulations of the individual 
spacing distribution in the bulk for the most relevant case (3=1 and various values of 
7, N, M and k (the location within the spectrum). Our findings are summarised in Fig. 
El displaying an excellent agreement between the numerical histogram and our surmise. 

In Fig. [6] (top), we keep 7 fixed to the value 7 = 7 and vary N,M and the 
location k within the spectrum. The histogram of the individual kth spacing reveals the 
independence of the spacing on the location k and on or M (sufficiently greater than 
2). In Fig. [6] (bottom), we keep A^, M fixed to the values (10, 15) respectively. Increasing 
7 (7 = 1, 7, 90) we can nicely see the convergence of both the numerical histogram and 
theoretical curves towards the WD surmise eq. P3l) . 
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6. Applications and Discussion 

Before comparing to actual data, let us put our results into the context of other models 
and highlight their features. 

From a RMT point of view, the first classification issue is whether or not the 
considered ensemble is invariant. Preserving invariance leads to a complete solvability 
for all spectral correlators in our approach as well as in [8]. On the other hand, non- 
invariant models as in [HI UHl [36] permit generically to compute only the macroscopic 
density, either implicitly or explicitly. A second issue is that of RMT universality classes, 
and we argued in section 13.11 that our macroscopic density is universal in a weak sense 
under (non- invariant) deformations using independent random variables. 

A third and important issue is to classify how the spectral density decays for large 
arguments, even though this may be difficult to appreciate from a small set of data. We 
have exemplified spectral densities with a power law and exponential decay using a x^- 
and inverse x^-distribution respectively. 

About the first issue above, we mention that the non-invariant ensemble [10] using 
a multivariate student distribution - a product of Gaussian and inverse x^-distributed 
random variables - leads to a density with power law decay. Another example is the 
deformed Gaussian Orthogonal Ensemble (GOE) [36] having an exponential tail as in 
our more general eq. (l32i) . for the special case of 7 = 2 (and c = 1). This model is 
non-invariant as well and can be derived from another entropy maximisation procedure 

The most natural setting to apply WL or its generalisations is in time series analysis. In 
Fig. [7] we compare to eigenvalues from financial covariance matrices for 2 different sets 
of data [3H1 139]^ . exploiting superstatistical models with x^- and inverse x^-distribution 
(displaying power law and exponential decay respectively). While in the top plot the 
former gives a better fit, in the lower plot our new eq. (l40l) for the inverse x^-distribution 
gives at least a comparable if not better fit to the data. We also compare to the standard 
MP distribution eq. fl38|) . which appears to give clearly a less good fit. Because in 
financial data there is no easy underlying physical principle for identifying extensive 
variables and thus the distribution class to be applied, our comparison must be heuristic, 
deciding case by case which distribution gives the best fit. 

A similar approach has been taken in a comparison with partly chaotic billiards 
|22j . Here the microscopic spacing distribution from various generalised WD classes is 
compared to the data, and the inverse x^-distribution gives the best fit. Other examples 
where the inverse x^-distribution has been recently identified is in turbulent fiows |41j . 
although not in a RMT setting. 

A very recent example of an application of a superstatistical RMT to complex 
networks was given in [37]. Here the deformed GOE model [36] within the class 
of exponentially decaying densities is compared to local statistics of the adjacency 

% Wc kindly thank the authors for permission to use the data from their papers. 
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Figure 7. A comparison to 2 sets of spectral densities from financial covariance 
matrices: top S&P500 price fluctuations of TV = 406 assets of M = 1309 days 
between 1991-1996 [40] (see also [38] Fig 1.), and bottom daily price exchanges of 
the Johannesburg Stock Exchange from Jan 1993-Dec 2002 [39]. In both cases a few 
higher eigenvalues are omitted. The smooth curves are the spectral densities for the 
given c < 1 from a ^^-distribution (green) [8], inverse x^-distribution (red) eq. (|40p 
and for comparison, the standard MP distribution eq. (|38|) (black). 

matrix. Due to its non-invariance the microscopic RMT predictions there are obtained 
from simulations. It would be very interesting to calculate these correlation functions 
analytically in the framework of our model with a more general exponential decay, but 
this is going beyond the scope of this article. 

7. Conclusions 

Using the recently proposed methods of superstatistics, we provided a systematic 
framework for generalising the Wishart-Laguerre ensembles of random covariance 
matrices, retaining the exact solvability of the original model. This is achieved by 
allowing the ensemble parameter - the inverse variance of the data matrix elements 
- to fluctuate from one sample to another according to a certain distribution /, and 
then averaging over /. We have given a new compact expression for the distribution 
of matrix elements in the particular case where the ensemble parameter has an inverse 
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X^-distribution. This distribution has appeared in the hterature when modelhng the 
volatihty of financial markets and thus it is of interest for the analysis of large arrays 
of data. Averaging over this distribution, we are able to express the spectral statistics 
of the generalised WL ensembles as an integral over the corresponding statistics for the 
standard WL ensemble. This is the key ingredient to solve the model exactly, for finite- 
N as well as in both the macroscopic and microscopic limits. In our model we have full 
control over the interplay between the deformation parameter 7, the matrix size and 
their respective asymptotic limits. 

We have computed exactly several spectral quantities, first deriving a generalised 
semi-circle and generalised Marcenko-Pastur density in the macroscopic large-iV limit 
for square and rectangular data matrices respectively: in both cases, we obtain an 
exponential tail. Secondly, after discussing in detail the level spacing distribution in 
WL ensembles and the applicability of the standard Wigner's surmise, we determined 
the microscopic level spacing distribution for all three j3 using a new, 7-dependent 
surmise, which exhibits a good agreement with numerical simulations. 

Our findings are illustrated via an application to financial covariance matrices where 
we make a comparison between fits resulting from a x^- and an inverse x^-distribution. 
It would be very interesting to find further applications, in particular to time series of 
other kinds of data. 

Acknowledgements: This work started while one of us (A.Y.A-M) was in a short visit 
to Brunei University. The hospitality of the members of Department of Mathematical 
Sciences is acknowledged. A. P. Masucci is thanked for providing a reference. Financial 
support by EPSRC grant EP/D031613/1, European Network ENRAGE MRTN-CT- 
2004-005616 (G.A.), and European Union Marie Curie Programme NET- ACE (P.V.) is 
gratefully acknowledged. We are very grateful to R. Kiihn for discussions and kindly 
sharing his data. 

Appendix A. Numerical Simulations 

Since it is not quite trivial to generate matrices with a non-standard Gaussian 
distribution we briefiy describe the algorithm in detail. We focused on matrices with 
real elements {/3 = 1) as an example, being the most relevant case for applications. 
The algorithm used is the following: 

(i) Draw a random variable ^ from the inverse-x^ distribution 

/(0 = r^-'e-^/Vr(7 + i); 

(ii) Generate a random M x N matrix X whose entries are normal variables with 
vanishing mean, and variance = 1/(27^); 

(iii) Generate the covariance matrix W = X^X; 
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(iv) Diagonalise W, obtaining its positive eigenvalues {n[ , . • • ^ I^n j stands for 
the £th sample); 

(v) Store: 

- in the matrix row Ve j the rescaled eigenvalues {^f^} (j = l,---,^) (where 

- in the matrix row Sej the bare spacings sf^ = fif^ — fifli {j = 2, . . . , N); 

(vi) Iterate the procedure R times, so that i = 1, . . . , R; 

(vii) Plot i) a normalised histogram of all the entries Vej (average spectral density) and 
ii) given a certain k between 2 and A^, a normalised histogram of all the entries 
of S in column k, normalised by the mean (1/-R) Xl^i '^f.fc (individual spacing 
distribution at location k). 

The plots of the normalised histograms are then compared with the theoretical results 
for finite or infinite N respectively (see Figs. [3lHl and [6|) . 
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